function [wp_decomp,wp_nc_decomp,wp_coll_decomp,wg_decomp,wg_nc_decomp,wg_coll_decomp] = DO_sim_decomp_func(dist_N,sim_decomp_size,param,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,opt_p_N,opt_p_O,opt_x_N,opt_x_O,opt_eta_U)

    %----- Specify Parameters -----%
    beta = param(2);            % Discount factor (factors in exit probability)

    d = param(4);               % Human capital depreciation rate
    lambda = param(6);          % Probability of search OTJ
    chi = param(7);             % Investment cost parameter

    delta =param(11);           % Separation rate
    h_sigma =param(12);         % Standard deviation of initial human capital distribution
    gamma=2;

    % Specify grids %
    hlb = 0.0;                 % Human capital grid lower bound 
    hub = 7.0;                 % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % ======================= Create discretized h distribution ==================== %
    h_mu = 1;                  % Mean of initial human capital distribution
    cdf_h = zeros(1,nh);       % CDF of initial h distribution
    pdf_h = zeros(1,nh);       % PDF of initial h distribution

    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end

    % Find index for next period human capital (h') if the agent only lets human capital depreciate
    ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
    for ind_h = 1:nh
        h = h_grid(ind_h);
        error_min = 100;
        for ind_hp = 1:nh
            hp = h_grid(ind_hp);       % hp is the value of h'
            error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
            if (error<error_min)
                error_min = error;
                ind_hp_min(ind_h) = ind_hp;
            end
        end
    end
    
    % Given dist_N - create cdf of h values and probability of being each phi type %
    pr_phi_1 = sum(dist_N(1,:))/sum(sum(dist_N(:,:)));
    
    cdf_h_N = zeros(2,nh);
    for ind_phi = 1:2
        for ind_h = 1:nh
            cdf_h_N(ind_phi,ind_h) = sum(dist_N(ind_phi,1:ind_h))/sum(dist_N(ind_phi,:));
        end
    end

    %% =============================================================================================================== %%
    %======================================== Simulation for Wage Moments and Results ==================================%
    %================================================================================================================= %%
    
    % Need to follow "respondents" from entrance because SS distribution does not save the eta values %
    nt = 5;               % Number of quarters to follow each individual - in this case we want to compare wage in quarter 1 with wage in quarter 5
    simul_d(sim_decomp_size,nt) = struct;
    rng(123456789);         % Set random number seed so results are exactly the same each time the code is ran   
    
    for id = 1:sim_decomp_size   % Loop over each individual in the sample
       
        for t = 1:nt        % Loop over each time period the individual is observed
            if (t==1) %%% First period observed, want to assign worker types based on the distribution of non-outsourced workers, then randomly assign half to actually be outsourced
                % Now assign phi type
                r1 = rand;
                if r1 < pr_phi_1
                    simul_d(id,t).phi_ind = 1;
                else
                    simul_d(id,t).phi_ind = 2;
                end
                ind_phi = simul_d(id,t).phi_ind;
                
                % Assign starting h value
                r1 = rand;
                simul_d(id,t).h_ind = 1;  % If r1 is not bigger than any CDF values, assign to the lowest slot
                for ind_h = 2:nh
                    if (r1>cdf_h_N(ind_phi,ind_h-1)) 
                        simul_d(id,t).h_ind = ind_h;
                    end 
                end
                ind_h = simul_d(id,t).h_ind;
                h = h_grid(ind_h);
                
                % Now assign job type randomly
                r1 = rand;
                if (r1<0.5)
                    simul_d(id,t).stat = 2; % outsourced
                else
                    simul_d(id,t).stat = 1; % non-outsourced
                end
                
                % Assume eta is value they would choose if searching with given characteristics %
                simul_d(id,t).eta = opt_eta_U(ind_phi,ind_h);
                eta = simul_d(id,t).eta;
                
                if (simul_d(id,t).stat == 1)
                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                    hp = h_grid(ind_hp);
                    simul_d(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                else
                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                    hp = h_grid(ind_hp);
     
                    simul_d(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
                end

            else % so t>1
                simul_d(id,t).phi_ind = simul_d(id,t-1).phi_ind;             % Individual's phi value is permanent
                ind_phi = simul_d(id,t).phi_ind;
                    
                if ((simul_d(id,t-1).stat)==1)        % If in N last period  
                    % Recall h from last period and apply any investment
                    ind_h = simul_d(id,t-1).h_ind;
                    simul_d(id,t).h_ind = opt_ind_hp_N(ind_phi,ind_h);
                    ind_h = simul_d(id,t).h_ind;
                    h = h_grid(ind_h);
                            
                    % N last period - no separation shock - does NOT search OTJ 
                    simul_d(id,t).stat = simul_d(id,t-1).stat;
                    simul_d(id,t).eta = simul_d(id,t-1).eta;
                    eta = simul_d(id,t).eta;
                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                    hp = h_grid(ind_hp);
                    simul_d(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                        - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));
                       
                else                                                  % If in O last period  
                    % Recall h from last period and apply any investment
                    ind_h = simul_d(id,t-1).h_ind; 
                    simul_d(id,t).h_ind = opt_ind_hp_O(ind_phi,ind_h);
                    ind_h = simul_d(id,t).h_ind;
                    h = h_grid(ind_h);

                    % O last period - no separation shock - does NOT successfully searchs OTJ
                    simul_d(id,t).stat = simul_d(id,t-1).stat;
                    simul_d(id,t).eta = simul_d(id,t-1).eta;
                    eta = simul_d(id,t).eta;
                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);  
                    hp = h_grid(ind_hp);
                    simul_d(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_h)) + U(ind_phi,ind_h) + (chi*(hp - h*(1-d)))^gamma...
                        - beta*delta*U_hat(ind_phi,ind_hp) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                        - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp));  
              
                end % End if for stat (when not dead or died before period start) last period
            end % End if t==1
        end % End loop over time periods in panel
    end % End loop over ids in simul_size
    
    %------------ Untargeted: Relative Wage Growth -----------------%
    count_growth = 0;
    count_growth_noncollege=0;
    count_growth_college=0;

    count_total = 0;
    count_wp_noncollege =0;
    count_wp_college = 0;
  
    for id = 1:sim_decomp_size
        count_total = count_total + 1;
        count_growth = count_growth + 1;
        if (simul_d(id,t).phi_ind==1)
            count_growth_noncollege = count_growth_noncollege + 1;
            count_wp_noncollege = count_wp_noncollege + 1;
        else
            count_growth_college = count_growth_college + 1;
            count_wp_college = count_wp_college + 1;
        end
    end % end loop over id
    
    X_wp = zeros(count_total,5);
    X_wp_noncollege = zeros(count_wp_noncollege,4);
    X_wp_college = zeros(count_wp_college,4);
    
    wage_vect = zeros(count_total,1);
    wage_vect_noncollege = zeros(count_wp_noncollege,1);
    wage_vect_college = zeros(count_wp_college,1);
        
    X_wg = zeros(count_growth,7);
    X_wg_noncollege = zeros(count_growth_noncollege,6);
    X_wg_college = zeros(count_growth_college,6);
    
    wg_vect = zeros(count_growth,1);
    wg_vect_noncollege = zeros(count_growth_noncollege,1);
    wg_vect_college = zeros(count_growth_college,1);
    
    count_total=0;
    count_wp_noncollege =0;
    count_wp_college = 0;
    count_growth=0;
    count_growth_noncollege=0;
    count_growth_college = 0;
    
    for id = 1:sim_decomp_size
        
        count_total=count_total+1;
        wage_vect(count_total,1) = log(simul_d(id,t).wage);
        if  (simul_d(id,t).stat==2)
            X_wp(count_total,1) = 1;
        end
        if  (simul_d(id,t).stat==1)
            X_wp(count_total,1) = 0;
        end
                
        if simul_d(id,t).phi_ind==2
            X_wp(count_total,2) = 1;                            % college
        else
            X_wp(count_total,2) = 0;                            % non-college
        end       
        X_wp(count_total,3) = h_grid(simul_d(id,t).h_ind);      % h
        X_wp(count_total,4) = log(h_grid(simul_d(id,t).h_ind)); % lnh
        X_wp(count_total,5) = 1; %constant
                
        % Save Results for Non-College Type
        if (simul_d(id,t).phi_ind==1)
            count_wp_noncollege = count_wp_noncollege + 1;
            wage_vect_noncollege(count_wp_noncollege,1) = log(simul_d(id,t).wage);
            % Indicator for outsourced or not
            if  (simul_d(id,t).stat==2)
                X_wp_noncollege(count_wp_noncollege,1) = 1;
            end
            if  (simul_d(id,t).stat==1)
                X_wp_noncollege(count_wp_noncollege,1) = 0;
            end
            % Controls
            X_wp_noncollege(count_wp_noncollege,2) = h_grid(simul_d(id,t).h_ind);      % h
            X_wp_noncollege(count_wp_noncollege,3) = log(h_grid(simul_d(id,t).h_ind)); % lnh
            X_wp_noncollege(count_wp_noncollege,4) = 1;                              % constant
        end
        % Save Results for College Type
        if (simul_d(id,t).phi_ind==2)
            count_wp_college = count_wp_college + 1;
            wage_vect_college(count_wp_college,1) = log(simul_d(id,t).wage);
            % Indicator for outsourced or not
            if  (simul_d(id,t).stat==2)
                X_wp_college(count_wp_college,1) = 1;
            end
            if  (simul_d(id,t).stat==1)
                X_wp_college(count_wp_college,1) = 0;
            end
            % Controls
            X_wp_college(count_wp_college,2) = h_grid(simul_d(id,t).h_ind);          % h
            X_wp_college(count_wp_college,3) = log(h_grid(simul_d(id,t).h_ind));     % lnh
            X_wp_college(count_wp_college,4) = 1;                                  % constant
        end

        count_growth = count_growth + 1;
                
        wg_vect(count_growth,1) = log(simul_d(id,5).wage) - log(simul_d(id,1).wage);
        % Indicator for outsourced or non-outsourced for both observations %
        if  (simul_d(id,1).stat==2 && simul_d(id,5).stat==2)
            X_wg(count_growth,1) = 1;
        end
        if  (simul_d(id,1).stat==1 && simul_d(id,5).stat==1)
            X_wg(count_growth,1) = 0;
        end

        % Indicator for college or non-college
        if simul_d(id,1).phi_ind==2
             X_wg(count_growth,2) = 1;                            % college
        else
             X_wg(count_growth,2) = 0;                            % non-college
        end
        % Other controls
        X_wg(count_growth,3) = h_grid(simul_d(id,5).h_ind);      % h
        X_wg(count_growth,4) = log(h_grid(simul_d(id,5).h_ind)); % lnh
        X_wg(count_growth,5) = simul_d(id,t).wage;
        X_wg(count_growth,6) = (simul_d(id,t).wage)^2;
        X_wg(count_growth,7) = 1;                               %constant
                
        if (simul_d(id,1).phi_ind==1) %----- Save data for non-college group only
            count_growth_noncollege = count_growth_noncollege + 1;
            wg_vect_noncollege(count_growth_noncollege,1) = log(simul_d(id,5).wage) - log(simul_d(id,1).wage);
            % Indicator for outsourced or non-outsourced for both observations %
            if  (simul_d(id,1).stat==2 && simul_d(id,5).stat==2)
                X_wg_noncollege(count_growth_noncollege,1) = 1;
            end
            if  (simul_d(id,1).stat==1 && simul_d(id,5).stat==1)
                X_wg_noncollege(count_growth_noncollege,1) = 0;
            end
            % Other controls
            X_wg_noncollege(count_growth_noncollege,2) = h_grid(simul_d(id,5).h_ind);      % h
            X_wg_noncollege(count_growth_noncollege,3) = log(h_grid(simul_d(id,5).h_ind)); % lnh
            X_wg_noncollege(count_growth_noncollege,4) = simul_d(id,1).wage;
            X_wg_noncollege(count_growth_noncollege,5) = (simul_d(id,1).wage)^2;
            X_wg_noncollege(count_growth_noncollege,6) = 1;                               %constant
                    
        else                        %----- Save data for college group only
            count_growth_college = count_growth_college + 1;
            wg_vect_college(count_growth_college,1) = log(simul_d(id,5).wage) - log(simul_d(id,1).wage);
            % Indicator for outsourced or non-outsourced for both observations %
            if  (simul_d(id,1).stat==2 && simul_d(id,5).stat==2)
                X_wg_college(count_growth_college,1) = 1;
            end
            if  (simul_d(id,1).stat==1 && simul_d(id,5).stat==1)
                X_wg_college(count_growth_college,1) = 0;
            end
            % Other controls
            X_wg_college(count_growth_college,2) = h_grid(simul_d(id,5).h_ind);      % h
            X_wg_college(count_growth_college,3) = log(h_grid(simul_d(id,5).h_ind)); % lnh
            X_wg_college(count_growth_college,4) = simul_d(id,1).wage;
            X_wg_college(count_growth_college,5) = (simul_d(id,1).wage)^2;
            X_wg_college(count_growth_college,6) = 1;    
        end
    end
    % Wage penalty regression -- No Sample Restrictions
    wp_decomp = regress(wage_vect,X_wp);
    
    % Wage penalty regression -- Only Non-College
    wp_nc_decomp = regress(wage_vect_noncollege,X_wp_noncollege);
    
    % Wage penalty regression -- Only College
    wp_coll_decomp = regress(wage_vect_college,X_wp_college);
    
    % Wage growth regression -- No Sample Restrictions
    wg_decomp = regress(wg_vect,X_wg);
    
    % Wage growth regression -- Only Non-College
    wg_nc_decomp = regress(wg_vect_noncollege,X_wg_noncollege);
    
    % Wage growth regression -- Only College
    wg_coll_decomp = regress(wg_vect_college,X_wg_college);

    

end % end function